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We present one-dimensional simulation results for the cold atom tunneling experiments by the 
Heidelberg group [G. Ziirn et al., Phys. Rev. Lett. 108, 075303 (2012) and G. Ziirn et al., Phys. 

Rev. Lett. Ill, 175302 (2013)] on one or two ®Li atoms confined by a potential that consists of an 
approximately harmonic optical trap plus a linear magnetic field gradient. At the non-interacting 
particle level, we find that the WKB (Wentzel-Kramers-Brillouin) approximation may not be used 
as a reliable tool to extract the trapping potential parameters from the experimentally measured 
tunneling data. We use our numerical calculations along with the experimental tunneling rates for 
the non-interacting system to reparameterize the trapping potential. The reparameterized trapping 
potentials serve as input for our simulations of two interacting particles. For two interacting (distin¬ 
guishable) atoms on the upper branch, we reproduce the experimentally measured tunneling rates, 
which vary over several orders of magnitude, fairly well. For infinitely strong interaction strength, 
we compare the time dynamics with that of two identical fermions and discuss the implications of 
fermionization on the dynamics. For two attractively-interacting atoms on the molecular branch, we 
find that single-particle tunneling dominates for weakly-attractive interactions while pair tunneling 
dominates for strongly-attractive interactions. Our first set of calculations yields qualitative but 
not quantitative agreement with the experimentally measured tunneling rates. We obtain quantita¬ 
tive agreement with the experimentally measured tunneling rates if we allow for a weakened radial 
confinement. 


I. INTRODUCTION 

Open quantum systems are at the heart of many physi¬ 
cal phenomena from nuclear physics to quantum informa¬ 
tion theory [I, 2]. In fact, all “real” quantum systems are, 
to some extent, open systems. Interactions with the envi¬ 
ronment cause decoherence, resulting in non-equilibrium 
dynamics. It is often simpler to design experiments that 
probe non-equilibrium physics than it is to design exper¬ 
iments that probe equilibrium physics. Conversely, the 
theoretical toolkit for describing systems in equilibrium 
is generally much farther developed than that for describ¬ 
ing systems in non-equilibrium. 

Ultracold atom systems provide a platform for real¬ 
izing clean and tunable quantum systems [3-6]. Over 
the past few years, much effort has gone into de¬ 
scribing non-equilibrium experiments that are accessi¬ 
ble, within approximate or exact frameworks, to the¬ 
ory. Notable experiments are the equilibration dynamics 
of one-dimensional Bose gases [7], the spin dynamics of 
dipolar molecules in optical lattices with low filling fac¬ 
tor [8], and the tunneling dynamics of effectively one¬ 
dimensional few-fermion systems [9, 10]. This paper fo¬ 
cuses on the latter set of experiments. Specifically, the 
goal of the present work is to describe the tunneling dy¬ 
namics of few-fermion systems, which are prepared in 
a well defined quasi-eigenstate (metastable state), into 
free space. We consider small systems and directly solve 
the time-dependent Schrodinger equation in coordinate 
space. As we will show, this approach provides a means 
to quantify the importance of the particle-particle inter¬ 
action, covering time scales from a fraction of the trap 
scale to thousands times the trap scale. Alternatively, 
one could adopt a quantum optics perspective and pur¬ 


sue a master equation approach. 

Tunneling is arguably the most quantum phenomenon 
there is: If the system was behaving classically, tunneling 
would be absent [11]. Tunneling plays an important role 
across physics, chemistry and technology. The scanning 
tunneling microscope [12], for example, nicely illustrates 
how a physics phenomenon, the tunneling of electrons, 
has been turned into a powerful practical tool (the imag¬ 
ing of materials). The a-decay, i.e., the decay of a "^He 
nucleus from a heavy nucleus, is an example discussed 
in most undergraduate physics texts (see, for example, 
Ref. [13]). The typical picture is to identify an effective 
reaction coordinate and to obtain the tunneling rate from 
a WKB analysis. While powerful, such treatments com¬ 
pletely neglect the effect of interactions. Interactions also 
play a crucial role in sorting out under which conditions 
electrons in light atoms tunnel sequentially or simulta¬ 
neously [14]. The two-particle system considered in this 
work has been realized experimentally and is the possibly 
simplest scenario that deals with a truly open quantum 
system (the atoms can escape to infinity) in which inter¬ 
actions (short-range atom-atom interactions) play a cru¬ 
cial role. As we will show, even for this relatively simple 
set-up, matching theory and experiment is a non-trivial 
task. Of course, two-particle tunneling has been investi¬ 
gated previously in this and related contexts [15-21]. 

The remainder of this paper is organized as follows. 
Section II introduces the Hamiltonian, the Heidelberg 
experiment and selected simulation details. Sections HI 
and IV discuss the molecular and upper branch tunneling 
dynamics. For both cases, it is argued that the trapping 
potential needs to be reparameterized. Using the repa¬ 
rameterized trapping potential, numerical simulations for 
the tunneling dynamics of two distinguishable ®Li atoms 
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on the molecular branch and the upper branch are dis¬ 
cussed. Comparisons with the experimentally measured 
tunneling rates are presented. Finally, Sec. V summa¬ 
rizes and provides an outlook. Simulation details and 
some technical aspects are relegated to Appendices A- 
G. 


II. SYSTEM HAMILTONIAN AND 
SIMULATION DETAILS 

A. One-body Hamiltonian, WKB analysis, and 
Heidelberg experiment 

This section considers a single ®Li atom with mass m. 
The atom is assumed to be in the hyperfine state |F', mj’). 
We consider the three lowest hyperfine states of the ®Li 
atom, referred to as |1) = |l/2,l/2), |2) = |l/2,—1/2), 
and |3) = |3/2,—3/2). Figure 1 shows the depen- 



a function of the magnetic field strength B. Solid, dashed and 
dotted lines correspond to states |1), |2), and |3), respectively 
(see text for details). States |1) and |2) are used in the upper 
branch experiments [9], while states |1) and |3) are used in the 
molecular branch experiments [10]. The higher-lying energy 
states shown by dash-dotted lines are not relevant for the 
present paper. 

dence of the hyperfine energy levels on the magnetic 
field strength B. The atom with coordinates {x,y,z) 
is trapped optically in a non-separable potential that is 
much tighter in the p-direction (p^ = x'^ + y'^) than in 
the z-direction [9, 10]. Throughout this work, we do not 
simulate the motion in the tight transverse confining di¬ 
rection. The transverse trapping frequency does, how¬ 
ever, enter into the calculation of the renormalized one¬ 
dimensional coupling constant (see Sec. IIC). Evaluat¬ 
ing the confinement created by the gaussian laser beam 
at p = 0, the effective one-dimensional single-particle 
Hamiltonian reads [9, 10] 

+ htrap(2) ■ZR, (i?)), (1) 

where the trapping potential Idrap along the z-direction 
depends implicitly on the internal or hyperfine state \j) 


TABLE I. Parameters from Refs. [9, 10] that define the trap¬ 
ping potential. Since the energy of the two-particle sys¬ 
tem on the molecular branch is smaller than the energy of 
the two-particle system on the upper branch, the p{t = 0) 
value for the molecular branch is chosen to be smaller than 
that for the upper branch; this guarantees that the tunnel¬ 
ing rates for the two experiments have roughly comparable 
orders of magnitude. The harmonic oscillator units are de¬ 
fined in terms oi uj = 2tt x 1234Hz, corresponding to Eho = 
8.177 X 10“®^ J, Oho = 1.167pm, and uj~^ = 1.290 x 10“'*s, or 
IJ = 1.223xlO®°Eho, Im = 8.570x 10®aho, and Is = 7753a;"h 
In an alternative levitation measurement, the magnetic field 
gradient was found to be B' — 1890(20)G/m ]9]. 


Quantity 

Value 

Go 

ks X 3.326pK = 56.16i?ho 

Zr 

9.975(5) x 10"®m = 8.548(5)aho 

p{-tr) 

0.795 

pit = 0) (upper branch) 

0.6875 

pit = 0) (molecular branch) 

0.63496 

dp/dt (for —tr < t < 0) 

-43s“^ 

B' (WKB approximation) 

1892G/m 


of the atom through the coefficient C\j), 

f/;rap(’2, t;_P, Z;r, (i?)) = 

p{t)Vo ( 1 - —- ) - pbC|j) (H)z. (2) 

The first term on the right hand side of Eq. (2) ac¬ 
counts for the optical confinement. Vq denotes the maxi¬ 
mum depth of the trap, p(t) a time-dependent parameter 
\p{t) < 1], and zr the Rayleigh range of the laser beam 
that produces the confinement. The second term on the 
right hand side of Eq. (2) is linear in z and makes the 
tunneling possible, pe is the Bohr magneton and C\j-^{B) 
depends on the hyperfine state, magnetic field strength 
and magnetic field gradient B', 

Cy^{B)=ci,^{B)B'. (3) 

Here, c\j){B) is a dimensionless parameter close to 1 (see 
below for details). Table I summarizes the trap parame¬ 
ters reported by the Heidelberg group [9, 10]; the parame¬ 
ters are obtained from a combination of measurement and 
WKB analysis. Figure 2 shows Vtrap for C = 1892G/m 
and three different values of p. The solid line shows 
the typical confinement at the beginning of the exper¬ 
iment while the dashed and dotted lines show typical 
confinements during the hold time of the upper branch 
and molecular branch experiments, respectively (see be¬ 
low for details). 

Since the trapping potential changes with time, there 
exists no set of units that characterizes the system 
equally well for all times. Throughout, following Ref. ]9], 
we choose w = 2 tt x 1234Hz to define the o scillator units 
Eho, Oho, and Tho: Eho = huj, Oho = \/^/{mio), and 
Tho = 2 ' k / uj . 

The confining potential Idrap has a local minimum at 
•Zinin and a local maximum at Zb- To gain insight into 
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FIG. 2. (Color online) The trapping potential, Eq. (2), for 
C = 1892G/m and three different values of p, p = 0.795 (solid 
line), p = 0.6875 (dashed line), andp = 0.63496 (dotted line). 
Vo and zb. are fixed at the values reported in Table I. 


the harmonic approximation, we expand V)rap around its 
local minimum and calculate the frequency ujimpip) of 
the harmonic term, 


^trap {p) 



p{t)Vo {z^ - 

(4 + 4i„)' 


(4) 


In the absence of the magnetic field gradient B', the 
minimum of Vtrap is located at Zmin = 0. For a finite 
magnetic field gradient, the local minimum Zmin depends 
on the parameters of the trapping potential. The fre¬ 
quency Wtrap(p) can differ notably from the frequency w 
and provides, in some cases, a more na tural unit. We 
define Fltrap(p) = ^trap(p), atrap(p) = \/fl/{mU}trap{p)), 
and rtrap(p) = 27r/a;trap(p)- Note that these units de¬ 
pend explicitly on p(t); correspondingly, we specify p(t) 
when we use these units. 

The single-particle tunneling dynamics is, to a good 
approximation, described by an exponential decay. 


.^sp,iii(0 — Fgp in(^ref) CXp[ Tsp(^ ^ref)]; (5) 


where T’sp,in(i) denotes the probability of finding the par¬ 
ticle inside the trap, the tunneling rate 7 sp is assumed to 
be constant, and tref is a reference time. Within the 
WKB approximation (see, e.g.. Ref. [22]), the tunneling 
rate reads 

7WKB ^ fWKB^^ (g) 


where the frequency tunneling coefficient 

T are given by 


and 


rWKB _ C Vtrap(.2;min,i— 0 ) 

2'Kh 


(7) 


T = exp ^-2 Vtrap(z)\dzj . (8) 

In Eqs. (7)-(8), Vtrap is the trapping potential with p{t = 
0) (see Fig. 3 for the time dependence of p), 2:niin.t=o is 


the z-value at which Ftrap with p(t = 0) takes its local 
minimum, and the WKB energy e of state n is found by 
the consistency condition 

J ^J2m[e - Vtrapiz)]dz = -|- 0 nh. (9) 

Here, Ze. 2 , and with z^p < Ze ,2 < -^e.s are 
the three solutions of e — Vtrap)^:) = 0 and n with 
n = 0,1,2,... denotes the order of the semiclassical 
“bound state” of the trap. In theory, one has Fsp.in(i) + 
T’sp,out(i) = Ij with the initial condition Psp,in(—^r) = 1- 
Here Psp,out(i) denotes the probability that the particle 
has left the trap. The inside and outside regions are de¬ 
fined through z < Zb and z > Zb, respectively, with Zb 
corresponding to the barrier position at time t = . 

We now briefly review the experimental sequence em¬ 
ployed by the Heidelberg group [9, 10]. The experiment 
prepared the atom in an “eigenstate” of the deep trap 
(p = 0.795 at < = —tr) and then lowered the barrier by 
decreasing p{t) over a time period tr- At time t = 0, p{t) 
reached its minimum. After a variable hold time fhoid, 
the barrier was ramped back up over a time period tj.. 
At time t = thoid + tr, the experiment monitored the 
fraction Psp,out(i) of the particle that had left the trap. 
To obtain F’sp,out(i), the experiment was repeated many 
times for each t = thoid + tr and the data were averaged 
[each individual experiment yields Hsp.out(thoid -f t^) = 0 
or Ij. The time sequence is sketched in Fig. 3. In the 



FIG. 3. (Golor online) Schematic of the time sequence of the 
experiment. After initialization of the system, the dimension¬ 
less parameter p(t) decreases from pit = —tr) to p{t — 0) 
with a rate dp/dt = —43 s~^, remains constant for thoid 
(thoid ^ tr), and increases to its initial value over the time pe¬ 
riod tr. The measurement is performed at the time thoid -ftr. 

experiment [10], the initial condition was F’sp,in(—^r) < 1 
due to non-unit state preparation fidelity. While this 
changes the overall normalization, it does not change the 
tunneling dynamics. 

The coefficients C|j^(H), and correspondingly the 
C\j)iB), depend on the magnetic field strength B, which 
is used to tune the atom-atom scattering length. The 
coefficients cy'^iB) can, at least in a first analysis, be 
obtained using the Breit-Rabi formula [23] (see Ap¬ 
pendix A). For state |3), the Breit-Rabi coefficient 
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is independent of the magnetic field. For states 
|1) and 12), the dependence of the Breit-Rabi coefficients 
on the magnetic field strength B is comparatively strong 
when B is small {B < 600G) and weak when i? —>■ oo 
{B > 600G). References [9, 10] did not use the Breit- 
Rabi formula to determine the c\j){B) coefficients (see 
below for details). 

To parameterize Vtrap, Refs. [9, 10] fed the result from 
“calibration measurements” into Eqs. (6) and (9). In a 
first step, the parameters Vq and zr of the optical trap, 
which is independent of the hyperfine state and magnetic 
field strength, were calibrated assuming p = 1. Specif¬ 
ically, the single-particle trap energy levels of the pure 
optical trap [C|j)(R) = 0 in Eq. (2)] were measured spec¬ 
troscopically and the parameters Vq and zr were chosen 
such that the WKB energy levels agreed with the mea¬ 
sured energies (see the supplemental material of Ref. [9]). 

For the upper branch tunneling experiment, p{t = —tr) 
and p{t = 0) were obtained by measuring the rela¬ 
tive integrated light intensities of the trap beams, i.e., 
p{t = —tr) and p{t = 0) were calibrated relative to 
p = 1 [24]. To obtain B', tunneling experiments at 
various magnetic fields using ®Li in state |2) were per¬ 
formed [25]. To prepare the atom in an excited trap 
state, the experiments used a trick. Two atoms in the 
same hyperfine states were prepared in the trap (these 
atoms do not interact), forcing the two-particle system 
to sit in a superposition of the lowest and first excited 
trap states. The assumption was then that the tunnel¬ 
ing dynamics proceeds as if there were a single particle 
in the first excited trap state and another single particle 
in the lowest trap state. The tunneling was attributed 
to the particle in the first excited trap state while the 
particle in the lowest trap state was assumed to have no 
chance of tunneling. This assumption is, as our simu¬ 
lations show, justified quite well (see Appendix B). To 
analyze the tunneling data, C|j) was assumed to be equal 
to 1 for all magnetic field strengths and B' was adjusted 
to yield a WKB tunneling rate 7^^^ that agreed with 
the measured tunneling rate The resulting B' was 

then used for all hyperfine states. 

The two-particle molecular branch experiments were 
conducted at magnetic field strengths varying from 350G 
to 1202G and utilized states |1) and |3) [10]. The param¬ 
eters p{t = —tr), Vb, zr, and B' were taken as those 
obtained from the upper branch experiments. Gompared 
to the upper branch experiments, p{t = 0) was reduced to 
obtain tunneling times smaller than a few thousand mil¬ 
liseconds and the magnetic field dependence of the coeffi¬ 
cients C|j) (B) was found to play a non-negligible role. For 
technical reasons, p{t = 0) was not calibrated via a “di¬ 
rect” photodetector measurement [24]. Instead, p{t = 0) 
and C|j) {B) were determined based on the WKB analysis 
of the experimentally measured single-particle tunneling 
rates (see supplemental material of Ref. [10]). Specif¬ 
ically, the single-particle tunneling measurements were 
performed at i? = 350G and 569G and the parameters 
p{t = 0) and c\j) (B) were adjusted to yield a WKB value 


7 ^^® that agreed with the measured tunneling rate 7 |pP 
at both 5-fields (see supplemental material of Ref. [10]). 
The analysis yielded p(f = 0) = 0.63496 [10]. The C|i)(il) 
and C| 3 )(i?) values are given in Table II. 

B. Simulation of single-particle tunneling dynamics 

To determine the single-particle tunneling rate theo¬ 
retically, we prepare the initial state {t < —tr) through 
imaginary time propagation. The initial state can be 
thought of as a quasi-eigenstate. We then propagate the 
initial state in real time for t > -tr- For —G < t < 0, 
we change p{t) according to dp/dt = —43s“^. For t > 0, 
p{t) is kept constant, i.e., p{t) = p(0). By analyzing the 
flux through z = Zf,, we calculate Psp,in(i) and Psp,out(i)- 
We do not simulate the up ramp, i.e., the time period 
^hoid < t < thoid + tr, since we found that the popula¬ 
tions F’sp,in(i) and Psp,out(i) do not change appreciably 
during the up-ramp. The simulation details are described 
in Appendices C and D. 

C. Two-body Hamiltonian and simulation of 
two-particle tunneling dynamics 

This section considers two ®Li atoms, each described by 
the single-particle Hamiltonian [see Eq. (1)], that in¬ 
teract through the short-range potential Vjnt(' 2 i 2 ), where 
Zi 2 = zi — Z 2 . The two-body Hamiltonian H reads 

H (zi, Z 2 , t;p, ZR, (H), C|j 2 ) (B)) = 

7J"P(zi,t;p,ZR,C|ji)(H)) -f 
7J"P(z2,t;p,ZR,C|j2)(B)) + VirA{zi2)- (10) 

Since the range of the true ®Li-®Li van der Waals po¬ 
tential is, for the experiments considered, much smaller 
than the de Broglie wavelength of the atoms, the de¬ 
tails of the interaction potential are not probed and the 
true interaction potential can be replaced by a simpler 
model potential that has the same three-dimensional s- 
wave scattering length osd as the true atom-atom poten¬ 
tial. Eor ®Li the most precise magnetic field dependence 
of osD is given in Ref. [26]. To convert osd to the one¬ 
dimensional coupling constant giD, we assume a three- 
dimensional zero-range potential and strictly harmonic 
confinement with angular frequency Wp in the tight di¬ 
rection. Describing the two-body interaction potential 
along the z-direction by 

Fzr(^12) = giT:id{zi2), (11) 

the renormalized one-dimensional coupling constant giD 
is given by [27] 

giD _ 2a3D / |C(l/2)|a3D \ ^ .^^2) 

SwpOp Op V y/2 ttp ) ’ 

where C(l/2) is equal to —1.46035 and Op denotes the har¬ 
monic oscillator length in the tight confining direction. 
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ap = ^Jh/{mujp). The one-dimensional coupling constant 
gm and the one-dimensional scattering length oid are re¬ 
lated via oid = —To determine Wp, Ref. [28] 
analyzed the optical single-particle trap with p{t) = 1 in 
the absence of the magnetic field gradient, accounting 
for the longitudinal (weak) and transverse (tight) direc¬ 
tions. The harmonic frequency ujp in the transverse di¬ 
rection was found to be = 27r x 14.22(35)kHz [28]. 

For p{t) ^ 1, needs to be multiplied by \/p{t), i.e., 
[9; 10; 28], resulting in a time-dependent 
giu- As discussed at the beginning of Sec. IIIB, the time 
dependence of (/id has a negligible affect on the tunnel¬ 
ing rate and we neglect it for the calculations presented 
in Secs. IIIB and IVB. 

The addition of the linear term [second term on the 
right hand side of Eq. (2)] moves the atoms away from 
the origin to positive z values. Using Eq. (3) of the sup¬ 
plemental material of Ref. [28] to model the confinement 
created by the gaussian beam in the longitudinal and 
transverse directions and expanding around p — 0, one 
finds that the harmonic frequency in the transverse direc¬ 
tion decreases with increasing z. For z = Zmin {z = Zb), 
we find that the harmonic frequency in the p-direction 
decreases by around 14% (38%) and 11% (44%) for the 
molecular and upper branches, respectively, compared 
to the frequencies for z = 0. This suggests that the 
tight confinement length Op in the presence of the mag¬ 
netic field gradient may be larger than [p{t = 0)]“^/"^ap®^, 

where Uref = ^/(ww™^), and correspondingly that the 

coupling constant (/id is modified. We return to this as¬ 
pect in Secs. IIIB and IVB. We reemphasize that the 
renormalization prescription given in Eq. (12) relies on 
the harmonicity of the confinement. It is well docu¬ 
mented in the literature that this renormalization pre¬ 
scription is modified by anharmonicities [29-31]. 

Eor the molecular branch, it has been shown theoret¬ 
ically that the strictly one-dimensional energies for the 
system without tunneling agree quite well with the full 
three-dimensional energies provided the one-dimensional 
scattering length oid is larger than the harmonic oscil¬ 
lator length Op [32]. Correspondingly, we restrict our 
molecular branch calculations to this regime (i.e., the 
smallest aip con sidered in Sec. IIIB—calculated using 
(jjp = \Jp{t = 0)a;p®*^—is oid = l-IISoho, corresponding 
to (/ID = —1.797aho)- It should be kept in mind, however, 
that the validity regime of the one-dimensional frame¬ 
work could be different for static (energies) and dynamic 
(tunneling) observables. 

We use two different model interaction potentials Vint, 
a zero-range potential Vzr, Eq. (11), and a finite-range 
gaussian potential Vfr, 

Vfr( 2 : 12 ) =-Vg exp , (13) 

where Vg and zq denote the depth (Vg > 0) and the 
range of the interaction. We use zq = 0.3aho, 0.2aho and 
O.lflho, and adjust Vg for each zq such that Vfr yields 


the desired one-dimensional two-body coupling constant 
(/ID- Throughout, Vg is chosen such that Vfr supports 
at most one even parity bound state in free space. We 
find that the dependence of the tunneling observables 
on the range zq is small. This together with the fact 
that Vzr and Vfr yield compatible tunneling results (as 
discussed below, we checked this for selected parameter 
combinations) justifies the use of comparatively large zq. 
The real time propagation of the two-particle system is 
discussed in Appendices C and E. 

To get a first feeling for the two-particle system, we 
consider the system with p = p{—tr) and map out the 
energy spectrum as a function of (/id- We use the 
imaginary time propagation (see Appendix D) to find 
the “eigenenergies” and “eigenfunctions” of the system 
[strictly speaking, the states are metastable due to the fi¬ 
nite barrier for p{—tr) = 0-795]- Eigure 4 shows the spec- 



FIG- 4. (Color online) Energies of two interacting trapped 
particles as a function of —1/giu- Solid and dashed lines show 
the energies for two particles with zero-range interaction and 
finite-range interaction, respectively, in a harmonic trap with 
frequency Wtrap- Circles and diamonds show the energies for 
two particles with zero-range interaction and finite-range in¬ 
teraction, respectively, in an anharmonic trap [see Eq. (2)]. 
Both particles feel the same external potential [p = 0.795, 
^\ 3 ){B) ~ 1.00115, and B' = 1892C/m, corresponding to 
(iJtrap = 27r X 1067.87Hz]. The width of the finite-range po¬ 
tential is 20 = 0.0930Otrap. 

trum for two interacting particles described by the Hamil¬ 
tonian H, Eq. (10), with p = 0.795 and zr = 8.548aho 
as a function of — 1/(/id- Both particles are assumed 
to feel the same single-particle trapping potential with 
C = 1894.18G/m. Diamonds and circles show the en¬ 
ergies for the zero-range potential and the finite-range 
potential with zq = O.loho = 0.0930atrap, respectively. 
Note, throughout we use the zero-range potential to de¬ 
scribe the positive (/id portion of the upper branch. In 
this regime, the Hamiltonian with finite-range interaction 
supports many deep-lying states, making it challenging 
to select the low-energy states of interest (recall, the rel¬ 
ative and center-of-mass degrees of freedom are coupled). 
Alternatively, one might consider using a purely repulsive 
finite-range two-body potential. In this case, however, a 
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large ^id would require a large range, thereby making 
the calculations mo del-dependent. Hence, this alterna¬ 
tive approach is not pursued here. Figure 4 uses the nat¬ 
ural units fltrap and i^trap with Wtrap = 27r X 1067.87Hz 
[see Eq. (4)]. The agreement between the zero-range and 
finite-range energies is very good for the gio considered. 

To illustrate the effect of the trap anharmonicity, solid 
and dashed lines show the eigenspectrum for two par¬ 
ticles interacting through Vzr and Hfr under external 
harmonic confinement with frequency wtrap (i-e., with¬ 
out magnetic field gradient and without anharmonic¬ 
ity) . The solid and dashed lines agree very well for most 
5 id- Differences are visible for the “diving” states near 
l/SiD ~ 0. The differences arise because the states with 
odd relative parity are not affected by the zero-range 
potential but are affected by the finite-range potential. 
Comparing the energy spectrum for the isotropic trap 
(lines) and the anharmonic trap (symbols), we see that 
the energies of the lowest state agree well for negative (/id 
( molecular branch) and positive gm (upper branch). The 
negative (/id portion of the upper branch is affected com¬ 
paratively strongly by the anharmonicity. In this regime, 
the anharmonicity leads to a decrease of the energies due 
to the widening of the trap. The coupling between the 
relative and center-of-mass degrees of freedom leads to 
avoided crossings between the energy levels that corre¬ 
spond, for the harmonic trap, to even relative and odd 
relative parity states. The eigenstates corresponding to 
the symbols on the upper branch and molecular branches 
serve as initial states for the real time evolution, i.e., 
these states serve as our initial wave packets at t = —tr- 

To analyze the tunneling dynamics of the two-particle 
system, we partition the configuration space as shown in 
Fig. 5. Region i ?2 corresponds to the situation where two 
particles are in the trap, region Rq corresponds to the sit¬ 
uation where both particles have left the trap, and region 
Ria {Rib) corresponds to the situation where particle 1 
( 2 ) has left the trap while particle 2 ( 1 ) is in the trap. 
The regions Riau, Ribu, and i?on correspond to numer¬ 
ical regions in which we apply damping (see below). The 
region Rj is encircled by the boundary Bj (the B^ ’s are 
not labeled in Fig. 5). To analyze the flux, the bound¬ 
aries Bj are broken up into boundary segments bjji that 
border regions Rj and Rj>. 

The flux through boundaries 62 ,lA and 62,15 is inter¬ 
preted as uncorrelated single-particle tunneling while the 
flux through boundary 62,0 is interpreted as pair tun¬ 
neling. The pair tunneling rate extracted from the flux 
through 62,0 is not unique and depends on Zpair- Sec¬ 
tion HIB considers l.loho < oid < 4.5aho; motivated by 
the fact that the size of the free space molecule is approx¬ 
imately oiD [33], we use 2 :pair = 2aiD for this oid range. 
For the upper branch simulations, we use 2 ;pair = 2aho- 
We found that the flux through 62,0 is vanishingly small 
for the upper branch simulations. We set Zi/o, which de¬ 
fines where the “inside” region ends and the “outside” 
region starts, such that > max(zt,i, 2 ^ 62 )■ For the 
molecular branch simulations, the flux dynamics is quite 


^ 2 ^ ^ 


p 

■'iBn 

O’ 

03 

3 

0 

3 


b ^On 

•^0,0n 

^IB 

^2,1B 

0 

03 

0 


Ro 

0 

p 

0 

3 


P 

0 

''v^pair 

^lAn.On 




^1A,0 


R2 



R 

5 ■'lA 

> M 

> ^ 


^bl ^i/o Zj 

FIG. 5. (Color online) Configuration space of the two-particle 
system. The regions Rj {j = 2, lA, IB,0, lAn, lBn,0n) are 
shown in different colors/shades. Each region Rj is sur¬ 
rounded by the boundary Bj (not shown). Boundary seg¬ 
ments that divide regions Rj and Rj/ are labeled by bjj/. Zbi 
and Zb 2 denote the position of the maximum of the barrier 
at t = —tr- Zi/o divides the “inside” from the “outside”; we 
choose Z\/o to be larger than Zbi and Zb 2 to ensure that the cal¬ 
culated flux is independent of how it is extracted. Zd denotes 
the largest z\ and 22 for which we calculate the “physical” 
wave packet. Zpair is equal to 2aiD for the molecular branch 
and equal to 2aho for the upper branch; Zpair enters into our 
analysis of the pair tunneling (see text for details). 

complex near the top of the barrier. To be independent 
of the “near-held” dynamics, we choose Zj/o ~ ISoho and 
13aho for the molecular branch and upper branch, respec¬ 
tively. The physical regions end at Zd, i.e., for zi > Zd 
or Z 2 > Zd a damping function is applied. The damp¬ 
ing function acts like an absorbing boundary (see Ap¬ 
pendix F). The damping function is needed since the 
fiux reaches the end of the simulation box within a small 
fraction of the total simulation time. Zd has to be so 
large that the two particles are essentially uncorrelated 
for z > Zd. In practice we vary Zd and choose its value 
such that the observables do not change as Zd is increased. 
Typical values for Zd are 25aho for the molecular branch 
simulations and 13aho for the upper branch simulations 
(for the upper branch, we found that Zd = zj/o yields 
the same results as Zd > Zj/o). As mentioned above, the 
time-dependent simulation starts at t = —tr, where the 
probability P 2 {—tr) to find two particles in the trap (i.e., 
in region R 2 ) equals I. For t > tr, ^ 2 ( 1 ) decays with time. 
This decay, except for a short period of time {t < 20ms), 
is well described by the exponential function 

P2{t) = P2{tref) exp[- 72 (t - Aef)], (14) 

where 72 denotes the decay rate. Since both uncor- 
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related single-particle tunneling and pair tunneling can 
contribute to the change of P 2 it), we break 72 into two 
pieces, 72 = 7 s + 7p, where 7 s and 7 p denote the single¬ 
particle tunneling and pair tunneling contributions, re¬ 
spectively (see Appendix G for details). A non-zero 7 s 
means that the probability Pi(t) to find one particle in 
the trap is finite. We also define the mean number N of 
trapped particles, 

N{t) = 2P2{t)+Pi{t). (15) 

The time dependence of N{t) is approximately parame¬ 
terized by an exponential decay with tunneling rate 7 , 

N{t) = fV(Gef) exp[- 7 (t - Gef)] + G, (16) 

where C denotes a constant. Subsections lllB and IV B 
present the results of our time-dependent two-particle 
simulations. 

III. MOLECULAR BRANCH TUNNELING 
DYNAMICS 

A. Single-particle tunneling dynamics and trap 
calibration 

In the following we perform exact numerical calcula¬ 
tions for the trap parameters reported in Table I. We 
will show that the numerically obtained tunneling rates 
do not agree with the measured ones and propose an al¬ 
ternative calibration approach. 

The trap employed in the molecular branch experi¬ 
ments was calibrated, in addition to the calibration ex¬ 
periments already discussed in Sec. II, based on four 
single-particle experiments [10] [see Table II and dia¬ 
monds in Figs. 6 (a) and 6 (b)]. In our first calculation, we 
use C|i) = I872.87G/m, corresponding to C|i) = 0.98989 
and B' = I892G/m, and prepare the system in the 
trap ground state [see the diamond in Fig. 6 (a)]. The 
dashed line in Fig. 7 shows the result of our simulation 
for p(0) = 0.63496. A fit of our data for t > I5ms (the 
short-time dynamics exhibits, as can be seen in the inset 
of Fig. 7, oscillations) to Eq. (5) yields 7 "™ = I5.39s“^ 
(see circles in Fig. 7). The tunneling rate 7 "™ ob¬ 
tained from the real time propagation is nearly twice 
as large as the experimentally measured tunneling rate 
Tip"", Tip"" = 8.28(0.49 )s-i [10]. This means that the 
trap parameters reported in Ref. [lO], obtained through 
the WKB analysis, yield a tunneling rate that deviates 
by a factor of nearly 2 from the experimentally measured 
tunneling rate, 7 ^p“/ 7 |pP = 1.86. To understand this, 
we treat as a parameter. Our tunneling simulations 
indicate that the exact shape of the initial state, and thus 
p{—tr)Vo, has a very small effect on the tunneling rate. 
The tunneling rate, in contrast, depends appreciably on 
the value of p(0)Vb- Thus, changing tr while keeping 
p{—tr)Vo fixed at 0.795Vb has a similar effect to changing 
Vb- Solid and dotted lines in Fig. 8 show the numerically 
determined tunneling rate and the WKB tunneling 



FIG. 6. (Color online) Parameter combinations {p{t — 
0),C|j)(i3)) that reproduce the experimentally measured 
single-particle tunneling rates at (a) B = 350G and (b) 569G. 
For all calculations, zn = 9.975/rm is used. In panels (a) and 
(b), the initial state corresponds to the trap ground state. 
The bands show the parameter combinations for which our 
full time-dependent calculations reproduce the experimentally 
measured tunneling rates. The widths of the bands origi¬ 
nate from the experimental error bars [10]. In panels (a) and 
(b), the dark (magenta) and light (cyan) bands correspond 
to ®Li atoms in states jl) and |3), respectively. Circles and 
squares show parameter combinations for states jl) and |3), 
respectively, that are used in the two-particle calculations (see 
Sec. IIIB). For comparison, the diamonds show the {p{t = 0), 
C|j)(R)) pairs that were suggested in Ref. [10]. 

rate 7 ^^® as a function of p{t = 0 ), i.e., for varying tr 
(using C|i) = 0.98989 and B' = 1892G/m [10]). It can 
be seen that the WKB analysis yields tunneling rates 
that differ by a factor of about 1 /2 from those obtained 
from the full time evolution. This is elaborated on fur¬ 
ther in Appendix H. Since the trap parameters reported 
by the experimental group utilized the WKB approxima¬ 
tion, we conclude that the trap parameters reported in 
Table I are inaccurate. Table II compares the measured 
tunneling rates 7 ™^ with the numerically calculated tun¬ 
neling rates 7 ”p ™ for the trap parameters summarized in 
Table I and the c\j) coefficients listed in Table II. 

The bands in Figs. 6 (a) and 6 (b) show the {p{t = 
0 ),C|j)) values for state | 1 ) [darker (magenta) band] and 
state 1 3) [lighter (cyan) band] for which the 7 ™™ agree 
with the experimentally measured single-particle tunnel¬ 
ing rates for states |I) and |3). In our calculations, the 
initial state corresponds to the lowest trap state. In a 
first attempt, we did set C|j)(i?) = c^^{B) and aimed to 
find unique values for p{t = 0) and B' that would repro¬ 
duce all four experimentally measured tunneling rates. 
For the functional form of the potential (with the pa- 





TABLE II. Experimentally measured single-particle tunneling rates for selected magnetic field strengths and initial single¬ 
particle states relevant to the molecular branch experiments [10]. Column 4 reports the values of the dimensionless coefficients 
cy){B) reported in Ref. [10]. The fifth column reports the tunneling rate 7 "p™ obtained from the exact time evolution using the 
trap parameters listed in Table I. As shown in column 6, the exact time evolution yields tunneling rates that are inconsistent 
with yjp^, suggesting that the trap calibration that involves the WKB analysis needs to be refined. 


state \j} 

BiG) 

(s"') 

c\piB) 

7“”“ (s“^) 

num /^exp 
Isp / /sp 

|1) trap’s gr. st. 

350 

8.28(0.49) 

0.98989 

15.39 

1.86 

|3) trap’s gr. st. 

350 

30.12(2.81) 

1.00311 

50.24 

1.67 

jl) trap’s gr. st. 

569 

21.76(1.12) 

0.99968 

37.36 

1.72 

|3) trap’s gr. st. 

569 

35.25(3.57) 

1.00457 

55.87 

1.58 



FIG. 7. (Color online) Single-particle tunneling as a function 
of time for a ®Li atom in state jl) ai B — 350G in the trap 
ground state. The dashed and solid lines show the probabil¬ 
ity Psp,in(t) of finding the particle in the trap calculated using 
the exact time evolution for p{t = 0) = 0.63496 and C|i) = 
1872.87G/m (c|i) = 0.98989 and B' = 1892G/m) (these pa¬ 
rameters are proposed in Ref. [10]) and for p{t = 0) = 0.63536 
and C|i) = 1862.03G/m (this is one of many parameter sets 
that reproduces the experimentally measured tunneling rate), 
respectively. The time evolution starts at —tr (tr ~ 3.72ms 
and tr ~ 3.71ms for the dashed and solid lines, respectively). 
Circles and squares show exponentially decaying functions 
[see Eq. (5)] with 7 sp = 15.39s“^ and 7 sp = 8.28s“^, respec¬ 
tively. The inset shows a blow-up of the short-time behavior. 


rameters Vb, zr, S', and dp/dt from Table I), such a 
parameter combination does not exist. Allowing zr to 
vary does not change the situation. To reproduce the ex¬ 
perimentally measured tunneling rates, we thus decided 
to treat Cy'j{B) as a free parameter. For example, we set 
C|3)(569G) = c^^(569G) and B' = 1890G/m and deter¬ 
mine pit = 0 ) such that we reproduce the experimental 
single-particle rate. Wefindp(t = 0) = 0.63536. We then 
set p{t = 0) to 0.63536 and find C|i)(569G), C|i)(350G), 
and C|3)(350G) such that 7 ™“ = 7 ™^ [see squares and 
circles in Figs. 6 (a) and 6 (b)]. We emphasize that these 
are not unique parameter combinations. Alternative pa¬ 
rameter combinations that are also used in Sec. Ill B are 
marked in Figs. 6 (a) and 6 (b). 

To obtain the Cy^iB) coefficients for other magnetic 
fields, we use interpolations/extrapolations. For state 



FIG. 8 . (Color online) Tunneling rate of a ®Li atom at 
B = 350G as a function of the dimensionless parameter 
p{t = 0). The atom is prepared in the ground state of 
the trap, and cp) = 0.98989 and B' = 1892G/m are used. 
The solid and dotted lines show the tunneling rates obtained 
through exact time propagation and the WKB approxima¬ 
tion, respectively. The horizontal band shows the tunneling 
rate 7 ^*^ = 8.28(0.49)s“^ measured experimentally [10] (the 
width of the band represents the experimental error bar). 

| 1 ), we use 

C|p(i3)«co + ^ + ^ (17) 

with Co = 1.00338, c_i = —1.89121G, and c _2 = 
— 1565.12G^. This functional form (i) reproduces 
C|i)(350G) = 0.985202 and C|i)(569G) = 0.995224 and 
(ii) is designed such that the functional dependence of 
c\yiB) is similar to that of Cp^(i3). For state |3), we use 
C| 3 )(^) = C|3)(569G) for B > 569G and a linear interpo¬ 
lation for 350G < B < 569G using the known C| 3 ) values 
at 350G and 569G. Table III summarizes the parame¬ 
ters that are used in Sec. Ill B to model the two-particle 
experiments. 

B. Two-particle tunneling dynamics 

This section considers two attractively-interacting ®Li 
atoms in hyperfine states |1) and |3) on the molecular 
branch. As discussed in Sec. IIC, the one-dimensional 
coupling constant giD depends on pit). Specifically, giD 
changes for t = —tr to t = 0 and is constant for t = 0 
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to t = ihoid- While this time dependence can be incorpo¬ 
rated straightforwardly into the finite-range simulations 
(in this case, the depth Vq can be made to vary with 
time), incorporating the time dependence into the zero- 
range calculations is more involved since giD enters into 
the propagator. To estimate the importance of the time 
dependence during the initial down ramp (time t = —tr 
to 0), we compared the simulation results for the cases 
where the full time dependence of giu w as a ccounted 
for [i.e., ujp was calculated according to and 

where the time dependence was neglected [i.e., Wp was 
calculated according to y^p(0)Wp®^] for selected magnetic 
field strengths. We found that the difference between 
the resulting tunneling rates is between 0.02% and 0.2%. 
Since this difference is much smaller than the difference 
between our calculated tunneling rates and the experi¬ 
mentally measured tunneling rates (see below), the time 
dependence of gio is neglected in what follows. The rea¬ 
son why the tunneling rates, calculated by accounting for 
and neglecting the time dependence of giu , are so similar 
is two-fold. First, very little tunneling occurs during the 
down ramp. Second, the overlap between the states at 
t = —tr with somewhat different giu is much larger than 
the overlap between the states at t = —tr and t = 0. This 
implies that the down ramp has a much larger effect on 
the state that results at t = 0 than a small variation of 
giD during the down ramp. 

The top panel in Fig. 9 shows the magnitude 
|j(zi,Z2,i)| of the flux for giu = —1.451iilhoaho [corre¬ 
sponding to ain/aho = 1-378 {B = 1202G, see Table III)] 
at t = 98ms. As can be seen (see also the arrows in the 
top panel of Fig. 9), the flux density is maximal along 
zi ~ 22- Only a small portion of the flux is directed along 
the Zi or Z2 directions. This demonstrates that pair tun¬ 
neling becomes dominant for sufficiently strong interac¬ 
tions. For comparison, the bottom panel of Fig. 9 shows 
the quantity Z 2 ,t)\ for the same t but giu = 0. In 
this case pair tunneling is absent. A careful compari¬ 
son of the flux in the Zi and Z2 directions shows that the 
flux along Z2 is notably larger, reflecting the fact that the 
trap felt by particle 2 (parameterized via C|3)) is shallower 
than the trap felt by particle 1 (parameterized via C|i))- 
We note that the flux has a very intricate structure in the 
vicinity of the barrier, especially in the upper panel, that 
is not visible on the scale of Fig. 9. Unlike the flux plots 
shown in Fig. 3 of Ref. [21], we do not observe “wave-like 
patterns” overlaying the flux. We speculate that these 
features are artifacts of the numerics of Ref. [21]. 

Figure 10 summarizes the tunneling rates obtained 
from our full time-dependent molecular branch simula¬ 
tions for finite-range interactions. To obtain these re¬ 
sults, ujp (and hence gio) was calculated according to 
— \/p(t = 0)Wp®^. Squares in Fig. 10(a) show the in¬ 
verse of 72, i.e., the inverse of the rate with which the 
probability P 2 {t) to find both particles in the trap decays, 
using the trap parameters that reproduce the experimen¬ 
tally measured single-particle tunneling rates. As can be 
seen, the squares lie notably above the experimentally 
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FIG. 9. (Color online) Probability flux \j{zi, Z2,t)\- The top 
and bottom panels show the probability flux at t = 98ms 
for two distinguishable particles with gin = —1.451 AhoUho 
and giD = 0, respectively (the trap parameters are given, 
respectively, in the sixth and first rows of Table III). The 
values of the flux are shown in the legend on the right in units 
of ij/aho (note the different scales for the top and bottom 
panels). The arrows indicate the primary directions of the 
flux j. 


measured (t™”^)”^ for finite giu] for gun = 0, the simu¬ 
lation results and the experimentally measured rate agree 
by construction since the single-particle tunneling rates 
in this case add up to 72 (see also Table III). 

The molecular branch tunneling dynamics has previ¬ 
ously been calculated by Landmark et al. using a time- 
independent method [21]. Unfortunately, the trap pa¬ 
rameters used to perform the calculations were not re¬ 
ported. The triangles in Fig. 10(a) show the result of 
this study. It can be seen that the inverse tunneling 
rate (7™™)”^ is a non-monotonic function of gun', such 
non-monotonic behavior is not displayed in our simula¬ 
tions. Reference [21] interpreted the non-monotonic de¬ 
pendence as an interplay between the trap parameters. 

To quantify the contribution of pair tunneling, we 
break 72 into two parts, 72 = 7p -b 7s, where 7p is 
the pair tunneling rate and 7s the single-particle tun¬ 
neling rate. We identify these rates from the flux pass¬ 
ing through the boundary 62,0 and the sum of the fluxes 
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TABLE III. Molecular branch dynamics for two distinguishable particles in states |1) and |3) for various magnetic field strengths. 
The second column reports the one-dimensional coupling constant giD calculated using LOp = •\/p(0)a;p'*^. The third column 
indicates whether the simulation results were obtained using the zero-range interaction model (ZR) or the gaussian interaction 
model with zo = 0.2aho (FR). Columns 4 and 5 report the Cy) coefficients for the trap parameterization with p{t = 0) = 0.63536 
(see Sec. II) and zb. = 8.548aho- Column 6 reports the tunneling rate 72 “™ [see Eq. (14)] obtained from our full time-dependent 
simulations. For comparison, column 7 shows the experimentally measured tunneling rates with error bars [34]. Column 8 
shows the rate 7 ^^ obtained from time-independent simulations ] 21 ]. 


B{G) 

yiD (<^ho-E^ho) 

ZR/FR 

C|i) (G/m) 

C| 3 > (G/m) 

7?™ (s-^) 

72 “" (s-i) [34] 

72 " (s-^) [21] 

569 

0 

— 

1881.11 

1891.32 

57.0 

57.01(3.74) 

— 

496 

-0.446 

ZR 

1877.16 

1890.19 

13.8(0.3) 

22 .2(1.0) 

19.2(0.5) 

496 

-0.446 

FR 

1877.16 

1890.19 

14.0 

22 .2(1.0) 

19.2(0.5) 

423 

-0.601 

FR 

1871.41 

1889.06 

6.67 

13.84(1.04) 

12.5(0.5) 

350 

-0.654 

FR 

1862.03 

1887.93 

4.27 

9.70(0.33) 

25.8(0.5) 

1202 

-1.451 

FR 

1891.38 

1891.32 

0.360 

2.14(0.19) 

0.4(0.5) 

1074 

-1.503 

FR 

1890.49 

1891.32 

0.293 

1.931(0.123) 

— 

958 

-1.595 

FR 

1889.44 

1891.32 

0.216 

1.227(0.053) 

— 

851 

-1.797 

FR 

1888.11 

1891.32 

0.137 

0.505(0.023) 

— 
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FIG. 10. (Color online) Molecular branch tunneling dynamics 
for two distinguishable particles as a function of gio- The cou¬ 
pling constant is calculated using utp = ^^/p(0)ujp‘^^. (a) The 
squares show the results from our full time-dependent sim¬ 
ulations using the trap parameters given in Table III; these 
trap parameters yield single-particle tunneling rates 7 "p ™ that 
agree with the experimentally measured single-particle tun¬ 
neling rates 7Jp*’. The symbols with error bars show the ex¬ 
perimental results [34]. For comparison, the triangles show 
the simulation results from Ref. ]21]. The inset shows the 
strongly-attractive region using the same symbols as in the 
main figure but a logarithmic j/-scale. (b) Squares show the 
ratio 7 p /72 obtained from our full time-dependent simula¬ 
tions. 


TABLE IV. Molecular branch dynamics for two distinguish¬ 
able particles in states ]1) and ]3) for various magnetic field 
strengths. The second column reports the one-di mensi onal 
coupling constant gm calculated using uip = 0.67^/p(0)ujp'‘^. 
The calculations are performed using the gaussian interaction 
model with zo = 0 . 2 aho and the trap parameters are the same 
as those for the calculations reported in Table III. Column 3 
reports the tunneling rate 72 *^™ [see Eq. (14)] obtained from 
our full time-dependent simulations. 


B (G) 

P'lD (^hoCf-ho) 

72™ («-') 

496 

-0.303 

22.4 

423 

-0.410 

13.4 

350 

-0.447 

9.57 

1202 

-1.018 

1.89 

1074 

-1.056 

1.53 

958 

-1.124 

1.04 

851 

-1.275 

0.62 


passing through the boundaries b 2 ,iA and 62 , 13 . Fig¬ 
ure 10 (b) shows the ratio 73/72 as a function of the in¬ 
teraction strength. We find that 73 is approximately 
equal to 0 for giu > —0.654i?hoaho- As one might pre¬ 
dict intuitively, the ratio 73/72 increases to close to 1 
for stronger attractive interactions. In this regime, the 
molecule can be treated as a point particle of mass 2 m. 
Our simulation results for giu > —0.654i?hoaho are con¬ 
sistent with the experimental observation of negligibe 
pair tunneling. In the strongly-interacting regime, i.e., 
for giD < —1.451i?hoaho, the experiments could not re¬ 
solve the pair versus single-particle tunneling fractions. 

To understand why our hnite-gio simulations predict 
larger tunneling constants I /72 than measured experi¬ 
mentally [see Fig. 10(a)], we repeated our simulations 
using several possible parameter sets that reproduce the 
experimentally measured single-particle tunneling rates, 
marked on the bands in Fig. 6 . We found that the 
two-body results remain almost unchanged, suggesting 
that the non-uniqueness of the trap parameterization is 
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not the cause for the disagreement. We also repeated 
one calculation using the zero-range interaction model 
as opposed to the finite-range interaction model (see Ta¬ 
ble III). Again, we found that the result remains almost 
unchanged, suggesting that finite-range effects are not 
the cause for the disagreement. As a third possibility 
we investigated the dependence of the tunneling rates on 
Wp. As we now show, a smaller ujp brings the tunneling 
rates obtained from the full time-dependent simulations 
in pretty good agreement with the experimentally mea¬ 
sured tunneling rates. 

As discussed in Sec. IIC the magnetic field gradient 
pushes the particles out to finite positive z, resulting in, 
on average, a weaker confinement along the tight con¬ 
finement direction. Squares in Fig. II show I/ 72 , ob¬ 
tained from our full time-dependent simulations, using 
the trap parameters that reproduce the experimentally 
measured single-particle tunne ling rate s and gijy calcu¬ 
lated according to Wp = = 0 ) 0 ;™*^ (see also Ta¬ 

ble IV). The factor of 0.67 yields (roughly) maximal 



SlD 

FIG. 11. (Color online) Molecular branch tunneling dynamics 
for two distinguishable particles as a function of gio. The cou¬ 
pling constant is calculated using ujp = 0.67^/p(0)uJp’‘^. The 
squares show the results from our full time-dependent simu¬ 
lations using the trap parameters given in Table III; these 
trap parameters yield single-particle tunneling rates 7"p™ 
that agree with the experimentally measured single-particle 
tunneling rates The symbols with error bars show 

the experimental results [34]. The inset shows the strongly- 
attractive region using the same symbols as in the main figure 
but a logarithmic t/-scale. 

agreement between the time constants obtained from our 
simulations and those measured experimentally (symbols 
with error bars in Fig. II). Recalling the discussion pre¬ 
sented in Sec. IIC, this value seems reasonable, though 
possibly slightly smaller than one might have expected 
naively. While other explanations for the disagreement 
between the squares and the symbols with error bars in 
Fig. 10(a) cannot be ruled out, our results indicate that 
the addition of the magnetic field gradient may have a 
non-trivial effect on the calculation of the renormalized 
one-dimensional coupling constant giu- 


IV. UPPER BRANCH TUNNELING 
DYNAMICS 

A. Trap calibration 

As discussed in Sec. II, the trap used in the upper 
branch experiment was calibrated by preparing two iden¬ 
tical non-interacting fermions in state | 2 ) at various mag¬ 
netic field strengths. The measured tunneling rates 7 ®’'p 
were obtained by fitting N{t) to an exponential plus 
a constant. Table V summarizes 7 ®’^p [35]. To see if 
the trap parameterization proposed by the experimen¬ 
tal group is accurate, we perform a time-dependent two- 
particle simulation for the anti-symmetrized two-particle 
wave packet using the trap parameters reported in Ta¬ 
ble I and C| 2 ) = 1. We find 7 ™™ = 6 . 86 s“^, which is 
about two times smaller than the experimentally mea¬ 
sured value, i.e., « 0.5 (note, this ratio is 

around 1.7 for the molecular branch; see Sec. Ill A and 
Appendix H). Similar to the molecular branch, we con¬ 
clude that the WKB approximation cannot be used to 
calibrate the trap. 

To recalibrate the trap, we set C| 2 ) = and adjust 
p{t = 0) and B' such that 7 "“™ for the anti-symmetric 
two-particle state at B = 782G agrees, within error bars, 
with the experimentally measured tunneling rate. As in 
the molecular branch (see Fig. 6 ), we do not find a unique 
parameter combination but a parameter band. Using 
p{t = 0) = 0.68, B' = 1890G and C| 2 ) = we find the 
tunneling rate 7 "“™ for several magnetic field strengths 
(see Table VI). Our 7 ™™ agree with 7 ®’'p within error 
bars, except for the cases at B = 750G and B = 855G, 
where the deviations are, respectively, about 1.1 and 2.5 
times larger than the error bars. 


B. Two-particle tunneling dynamics 

This section discusses the upper branch tunneling dy¬ 
namics for two distinguishable particles with finite inter¬ 
action strength giu- Solid and dashed lines in Fig. 12(a) 
show the mean number of trapped particles TV, Eq. (15), 
extracted from our full time-dependent simulations as a 
function of the hold time for two distinguishable parti¬ 
cles at B = 782G (gip = 192ahoAho; in what follows, 
we use giD = 00 for this magnetic field strength) and 
B = 900G {giu = —3.15ahoAho)- Here, gin is calculated 
using ujp = \/p{t = 0 ) 0 ;))®^. As can be seen in Fig. 4, the 
upper branch energy of the quasi-eigenstate at t = —U 
is larger for negative gm than for infinitely large ^id- 
This implies that the effective barrier height that the 
two-particle system sees is smaller at i? = 900G than at 
B = 782G, resulting in faster tunneling dynamics for the 
system at B = 900G than at B = 782G. The tunneling 
rates 7 , obtained by fitting our data to Eq. (16) or from 
the flux analysis (see Appendix G), are 7 ™™ = 127s“^ 
for B = 900G and 7 "^^“ = lbs"! for B = 782G. These 
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TABLE V. Upper branch dynamics for two distinguishable particles in states |1) and |2) for various magnetic field strengths. 
The second column reports the one-dimensional coupling constant gio calculated using tOp = Y^p(0)cJp®^. The third column 
indicates whether the simulation results were obtained using the zero-range interaction model (ZR) or the gaussian interaction 
model with zo — 0.2aho (FR)- Columns 4 and 5 report the C\j) coefficients; as discussed in Sec. IV A, we use c\j) = 

B' = 1890G/m, p{t = 0) = 0.68, and zb. = 8.548aho. Column 6 reports the tunneling rate 7“™ [see Eq. (16)] obtained from 
our full time-dependent simulations. For comparison, column 7 shows the experimentally measured tunneling rates with error 
bars [35]. 


B{G) 

plD (ciho^ho) 

ZR/FR 

C|i) (G/m) 

C| 2 ) (G/m) 

ynum (g-l) 

7 “" (s-') 

750 

6.15 

ZR 

1883.86 

1881.60 

4.2(0.5) 

2.9(0.2) 

782 

CX3 

ZR 

1884.56 

1882.47 

15 

12.8(1.2) 

855 

-4.42 

FR 

1885.88 

1884.07 

77 

62.8(8.2) 

900 

-3.15 

ZR 

1886.57 

1884.90 

127 

107(12) 

900 

-3.15 

FR 

1886.57 

1884.90 

130 

107(12) 


TABLE VI. Tunneling dynamics for two identical particles 
in state |2) for various magnetic field strengths. The second 
column reports the C\j) coefficients; as discussed in Sec. IV A, 
we use C|j) = B' = 1890G/m, p(t = 0) = 0.68, and 
Zb = 8.548aho. Column 3 reports the tunneling rate 7 "“™ [see 
Eq. (16)] obtained from our full time-dependent simulations. 
For comparison, column 4 shows the experimentally measured 
tunneling rate 7 ®^^ with error bars [35]. 


B{G) 

C| 2 ) (G/m) 

ynum (g-1) 

yexp (g-1) 

750 

1881.60 

13.2 

14.7(1.3) 

782 

1882.47 

13.8 

13.2(1.1) 

820 

1883.37 

14.5 

13.1(1.4) 

855 

1884.07 

15.1 

11.5(1.3) 

900 

1884.90 

15.8 

16.0(1.1) 


tunneling rates agree at the two sigma level with the ex¬ 
perimentally measured rates of = 107(12)s“^ [see 
triangles in Fig. 12(a)] and 12.8(1.2)s“^ [35] [see squares 
in Fig. 12(a)]. 

An important aspect of the tunneling dynamics of the 
upper branch is that the mean number of trapped parti¬ 
cles N decreases from 2 to approximately 1 over the hold 
times considered. This suggests that the particle that 
remains trapped has such a small energy that its tun¬ 
neling dynamics is orders of magnitude slower than the 
tunneling dynamics considered in Fig. 12. Indeed, we 
observe essentially no flux through the boundaries biA,o 
and biB,o- Comparing the portion of the wave packet in 
region Ria (or Rib) with the quasi-eigenstate of a single 
trapped particle shows that the remaining particle occu¬ 
pies to a good approximation the lowest trap state. This 
implies that the particle that leaves the trap carries away 
the “excess energy”. Performing single-particle calcula¬ 
tions for particles jl) and |2) initially in the trap ground 
state, we find tunneling rates of 0.008s“^ and 0.007s“^. 
This confirms the separation of time scales alluded to 
above. 

Circles in Fig. 12(b) show our tunneling time constants 
(ynu™)-! for two distinguishable particles as a function 
of the magnetic field strength. Our (7™™)“! follow the 
overall trend of the experimentally measured (y^’^P)"! 


[diamonds in Fig. 12(b)] but lie a bit lower (see also 
Table V). The discrepancy is largest for positive gin 
[B = 750G), where the dynamics is slowest. This is 
the regime where our simulations are, due to the slow 
tunneling, the most demanding. We estimate, however, 
that our numerical uncertainties do not account for the 
45% discrepancy between the calculated tunneling con¬ 
stant (7''“™)“^ and the experimentally measured tunnel¬ 
ing constant (7 ®’^p)“^. 

Motivated by the analysis presented in Sec. IIIB, one 
may ask how the tunneling rates for the upper branch 
depend on the Wp value used to calulate giu- We esti¬ 
mate that a scaling factor of around 0.85 improves the 
agreement between our simulations and the experiment 
for B = 750G; at the same time, the agreement for 
B = 855G and B = 900G detoriates. The fact that 
the “optimal” scaling factor for the upper branch seems 
to differ from that for the molecular branch is not unrea¬ 
sonable. First, since the non-linear trap term is larger, 
one might expect that u>p is modified less by the mag¬ 
netic field gradient term for the upper branch than for 
the molecular branch. Second, the excited upper branch 
states may be affected differently than the molecular 
branch states [one should keep in mind that Eq. (12) 
is an approximation]. 

It is interesting to compare, as has been done in 
the experiments, the tunneling dynamics for two distin¬ 
guishable particles with that for two identical particles, 
since two distinguishable but otherwise identical particles 
with infinitely large guy are known to become fermion- 
ized [27, 36, 37]. In the present case, the distinguishable 
particles in states jl) and |2) feel slightly different trap¬ 
ping potentials. Thus the fermionization concept does, 
strictly speaking, not apply. However, since C\i) and C\ 2 ) 
at B = 782G differ by only 0.2%, a meaningful compar¬ 
ison can be made. The dotted line in Fig. 12(a) shows 
the mean number of particles for two identical fermions 
in state |2). Since Cj2)(782G) < C|i)(782G), implying a 
higher barrier for the atom in state |2) than the atom 
in state jl), the non-interacting identical fermion sys¬ 
tem (two atoms in state |2)) tunnels slightly slower than 
the two distinguishable atom system (one atom in state 
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a “softer” confinement than the particle in state | 2 ), i.e., 
Wtrap for state |1) is less than Wtrap for state |2). Impor¬ 
tantly, the particles in states |1) and |2) at B = 782G 
( 51 D = 00 ) cannot pass through each other. Thus, since 
the particle in state | 1 ) tunnels slightly faster than the 
particle in state | 2 ) (see below) the probability Pzi>z 2 to 
have zi > Z 2 inside the trap gets depleted faster than 
the probability to have zi < Z 2 - Indeed, at t = 94ms, 
we have Pz^<z 2 = Pzi>z 2 - At the end of the simulation 
{t = 350ms), the probabilities to find an atom in state 
|1) and an atom in state |2) inside the trap are 48% and 
52%, respectively. 

In the “ideal fermionization scenario”, in which the 
infinitely strongly interacting particles feel the same ex¬ 
ternal potential, the ground state is two-fold degenerate. 
In our case, this degeneracy is broken since Cp) 7 ^ Cp). 
Solid and dotted lines in Fig. 13(a) show |'I'rei( 2 i 2 )P, 

pOO 

'^relizi2)= / 'i>{zi,Z2,t = -tr)dZcM, (18) 


FIG. 12. (Color online) Upper branch tunneling dynamics, 
(a) The dashed and solid lines show the mean nnmber of 
trapped particles N obtained from our full time-dependent 
simulations as a function of time for two distinguishable par¬ 
ticles at i? = 900G and B = 782G, respectively, using the 
trap and interaction parameters given in Table V. For com¬ 
parison, triangles and squares with error bars show the cor¬ 
responding experimental results [35]. The dotted line shows 
the mean number of trapped particles N obtained from our 
full time-dependent simulations as a function of time for two 
identical fermions at B — 782G, using the trap parameters 
given in Table VI. For comparison, circles show the corre¬ 
sponding experimental results, (b) Circles and triangles show 
the time constant 7 “^ obtained from our full time-dependent 
simulations for two distinguishable particles and two identi¬ 
cal fermions, respectively, as a function of the magnetic field 
strength B. For comparison, diamonds and squares with error 
bars show the corresponding experimental results [35]. 


where Zqm = (zi + Z 2 )/ 2 , for the ground state and the 
first excited state, respectively. The difference of the am¬ 
plitudes for Z 12 < 0 and Z 12 > 0 reflects the asymmetry 
of the trap potentials (see discussion above). The ground 
state wave function is greater or equal to zero everywhere 
while the first excited state wave function changes sign at 
Z 12 = 0. The energy difference between the two states is 
approximately 3 x lO^'^F'ho, corresponding to a time scale 
of 430ms. Since parity is not a conserved quantity and 
since the relative and center-of-mass degrees of freedom 
couple, we expect oscillations between the ground state 
and the first excited state at this time scale. Figure 13(b) 
shows the normalized overlap C>n„incM j 


O 


raremCM 


(t) 


('i'(Zi,Z2,t)l(/> 

TtrelT^CM (Z 1 ,Z 2 )} 


, ( 19 ) 


jl) and one atom in state | 2 )) with infinitely large ^id- 
Triangles and squares in Fig. 12(b) show the tunneling 
constants 7 “^ for two identical fermions as a function of 
P obtained from our simulations (see Sec. IV A and Ta¬ 
ble VI) and from experiment, respectively. Although the 
fermionization is only approximate. Fig. 12(b) shows that 
the tunneling rate curves for two distinguishable parti¬ 
cles and two identical non-interacting fermions cross at 
approximately B = 782G, corresponding to ^id = 00 for 
the | 1 )-| 2 ) interaction. 

Another consequence of the fact that C|2)(782G) < 
C|i)(782G) is that the probability to find the particle or¬ 
dering zi < Z 2 (or Zi > Z 2 ) for two atoms in states jl) 
and | 2 ) changes as a function of time. At t = —the 
probability Pzi>z 2 to find zi > Z 2 is 0.525 and the prob¬ 
ability Pzi<z 2 to find zi < Z 2 is 0.475 [see Fig. 13(a)]. 
This is due to the fact that the particle in state jl) feels 


between the time-evolving wave packet 4 '(zi,Z 2 ,t) 
and the two-body harmonic oscillator eigenstates 
- 2 ^ 2 ) with trap frequency wtrap and relative 
and center-of-mass quantum numbers n^ei and ncM- The 
solid line shows the overlap for the anti-symmetric refer¬ 
ence wave function (jjn.^incM with (nrepncM) = ( 1 , 0 ), 
which has odd relative parity. Dotted lines show the 
overlaps for states with even relative parity (see figure 
caption). The oscillation period, T r; 270ms, is compa¬ 
rable to but smaller than the estimated value of 430ms 
because the system is modified after t = —tr- Figure 13 
demonstrates that two distinguishable particles with in¬ 
finite 5 id but C|i) 7 ^ C| 2 ) exhibit unique dynamics that 
is absent for two identical fermions. It could be interest¬ 
ing in future work to tune the system toward and away 
from the ideal fermionization regime and to explore the 
resulting dynamics. 
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t (ms) 


FIG. 13. (Color online) Analysis of the upper branch time 
dynamics for two distinguishable particles interacting through 
a zero-range potential with gio ~ oo (the trap parameters 
are given in row 2 of Table V). (a) The solid and dotted lines 
show the density |4'rei(2i2)|^ of the lowest “eigenstate” and 
the first excited “eigenstate” at t = —tr. These states are 
nearly degenerate, (b) The dotted lines show the normalized 
overlap dnremcM between the wave packet 4'(zi, Z 2 ,t) and the 
harmonic oscillator states with even relative parity [from top 
to bottom, (nrei,n-cM) = (0,0), (2,0), (0,1), and (0,2)]. The 
solid line shows the overlap between the wave packet 

Z 2 ,t) and the lowest harmonic oscillator state with odd 
relative parity, i.e., with (nrei,ncM) = (1,0). The harmonic 
oscillator states are characterized by atrap = 1.073aho. 

V. SUMMARY AND OUTLOOK 

This paper provided a detailed analysis of the two- 
particle tunneling dynamics out of an effectively one¬ 
dimensional trap. Our studies were motivated by ex¬ 
periments by the Heidelberg group and our analysis was 
based on full time-dependent simulations of single- and 
two-particle systems. We found that the trap calibra¬ 
tion via a WKB analysis leads to an inaccurate trap 
parameterization; this finding is in agreement with a 
study by Lundmark et al. [21]. Using the reparameterized 
trapping potential, our tunneling rates for two identical 
fermions agree with the experimental results for all but 
two magnetic field strengths considered. 

Our simulations for the interacting two-particle sys¬ 
tems made a number of simplifying assumptions. The dy¬ 
namics in the tight confinement direction was only incor¬ 
porated indirectly via the renormalized one-dimensional 
coupling constant. For this, a harmonic trap in the tight 


direction was assumed. Moreover, we assumed simple 
short-range or zero-range interaction potentials. Deep- 
lying bound states and coupled channel effects were ne¬ 
glected entirely. Using the renormalized one-dimensional 
co uplin g constant giu with the transverse frequency 
^yp{0)uJp^^ as input, our simulations reproduced the up¬ 
per branch tunneling dynamics of the interacting two- 
particle system reasonably well. Our simulation results 
for the molecular branch dynamics agreed with the over¬ 
all trend of the experiment but did not yield quantitative 
agreement. We argued that the actual transverse confine¬ 
ment felt by the atoms in the presence of the magnetic 
field gradient may be weaker than in the absence of the 
magnetic field gradient. This motivated us to calculate 
the one-dimensional coupling constant using a weaker 
transverse trapping frequency as input. The resulting 
two-particle tunneling rates are in agreement with the 
experimentally measured rates over the entire range of 
magnetic field strengths considered. We note that our 
finding is consistent with Ref. [38], which found that the 
non-separability of a gaussian trap affects the tunneling 
rate in a double-well geometry. 

Our work suggests a number of follow-up studies. It 
would be interesting to extend the dynamical simulations 
to more particles and/or to include the tight confining 
directions. It would also be interesting to prepare other 
initial one- and two-particle states. For example, it would 
be interesting to investigate the tunneling dynamics from 
initial excited metastable states. 
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Appendix A: State dependence of the trapping 
potential and Breit-Rabi formula 

We consider an atom with total (orbital and spin) elec¬ 
tronic angular momentum quantum number J = 1/2 and 
nuclear spin I {I = 1 for ®Li). In the absense of an exter¬ 
nal magnetic field, the energy difference AW between the 
hyperfine states |F = J—1/2, mi?) and |F = 7-|-l/2,mF) 
is independent of mp- For ®Li with \F = 1/2) and 
|F = 3/2), AW is equal to 228.205 MHz [39]. According 
to the Breit-Rabi formula [23, 40], the energy {B) 

of the hyperfine state jU, tof) in an external magnetic 
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field of strength B is 


= - 


AVF 


±- 


AVF 


2(27+1) 

'imp 

^ 27 + 1 ^ 


giUBmpB 
1/2 


(Al) 


where x = {gj — gi)gBB/AW, gj is the Lande factor, 
and gi characterizes the magnetic moment of the nucleus. 
The plus and minus signs refer to states F = 7 + 1/2 
and F = 7—1/2, respectively. The constants gj = 
2.0023019(24) and gj = -0.0004476493(45) are deter¬ 
mined experimentally [41]. Figure 1 shows the magnetic 
field dependence of the hyperfine states of ®Li for F = 1/2 
and F = 3/2. The slope of these energy curves equals 
the negative of the magnetic moment of the atom [40], 
yielding 

= (A2) 

Equation (A2) characterizes the state and magnetic field 
dependence of the trapping potential (see Sec. II of the 
main text). The coefficients calculated according to 
Eqs. (Al) and (A2) are referred to as Breit-Rabi coef¬ 
ficients in the main text. 


Appendix B: Time dynamics for two identical 
fermions in an anti-symmetric state and in a 
product state 

The wave packet of two identical fermions is anti¬ 
symmetric under the exchange of the particles. To cali¬ 
brate the trap (see Sec. IV A), the assumption in using 
the WKB approximation was that the dynamics could 
be described as if a single particle was tunneling out of 
the first excited trap state. Our numerical simulations 
show that the tunneling rates are, indeed, very simi¬ 
lar. For the parameters listed in the third row of Ta¬ 
ble VI, we find 7 = 13.8s“^ for the two-particle system 
and 7 sp = 13.5s“^ for the single-particle system. As we 
discuss now, the tunneling dynamics is, however, quite 
different. 

Figure 14(a) shows the normalized overlaps C7n„incM 
[see Eq. (19)] between the time-evolving anti-symmetric 
two-particle wave packet Z 2 ,t) and the two-body 
harmonic oscillator eigenstates (/>nreincM(' 2 i) ^ 2 ) with trap 
frequency Wtrap and relative and center-of-mass quantum 
numbers rirei and ncM- The solid, dashed and dotted 
lines show the overlaps for n^ei = 1 and ncM = 0 , 1 , and 
2, respectively. The normalized overlaps oscillate for a 
short time {t < 10 ms) and quickly approach constants. 
We see essentially constant overlaps till the end of our 
simulation at f = 500ms. This indicates that the shape 
of the wave packet in region R 2 is constant in time. The 
overlaps vanish for even rirei, indicating that the anti¬ 
symmetry of the wave packet is preserved during the time 
evolution. 



t (ms) 

FIG. 14. (Color online) (a) Analysis of the time dynamics 
for two identical fermions. The solid, dashed, and dotted 
lines show the normalized overlap On^incM between the wave 
packet 4/(^1, Z 2 ,t) and the harmonic oscillator states with odd 
relative parity [from top to bottom, (nrei,ncM) = (1,0), (1,1), 
and (1,2)]. (b) Analysis of the time dynamics for a single atom 
in state |2), prepared in the first excited trap state at t = —tr. 
The lines show the normalized overlaps On{t) between the 
wave packet 4 ^( 2 , t) and the harmonic oscillator states with 
n = 0 — 4. For both panels, the parameters p{t = 0) = 0.68, 
2 r = 8.548aho, c™ = 0.99601, and B' = 1890G/m are used. 


Figure 14(b) shows the normalized overlap On(t), 


v'(4'(^,f)|T(z,f)) 


(Bl) 


between the time-dependent single-particle wave packet 
4f(z,f) and the time-independent single-particle har¬ 
monic oscillator functions (/)„(z) with quantum number 
n, n = 0 — 4. The overlaps shown in Fig. 14(b) oscillate 
at a frequency that is close to the natural trap frequency 
wtrap for t > 20ms. Moreover, the “envelopes” of the 
overlaps change in time, indicating that the shape of the 
wave packet in region R 2 changes with time. At t = 0, 
the wave packet has a finite overlap with the ground state 
harmonic oscillator wave function due to the change of 
the trapping potential. The contribution of the harmonic 
oscillator ground state to the wave packet is almost con¬ 
stant in time while the contributions of higher energy 
states deplete. This results in the increase of the nor¬ 
malized overlap oo{t) [see Fig. 14(b)]. In other words, as 
Tsp.in(f) decreases, the wave packet looks more like the 
lowest-lying trap state as opposed to the initial state. As 
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a consequence, the decay of Psp,in(0 with time deviates 
slightly from an exponential. 


Appendix C: Time propagation via Chebyshev 
expansion 

The time evolution of the two-particle wave packet 
is given by 

'll{zi,Z2,t) =U{t - to)'^{zl, Z2,to), (Cl) 

where the time-evolution operator — tg) is 

U{t — to) = exp[—iH{t — to)/h]. (C2) 

To evaluate Eq. (Cl), one has to expand the time- 
evolution operator U{t — to) in powers of —iH{t — to)/h. 
It has been shown that expanding U{t — to) in terms of 
the complex Chebyshev polynomials (j)k, 

Hit -to) = J2 ~ ) , (C3) 

fc =0 ' 

provides an efficient means to determine the time evo¬ 
lution of the wave packet [42]. Here, i? is a real and 
positive number that has been introduced to normalize 
the argument of (fik such that —iHit — to)/ihR) G [— I,*]- 
A key point is that the recursion relation 

(t>kiX) = 2X4>k-i{X) + 4>k-2iX) (C4) 

for the fc**' Chebyshev polynomial enables one to readily 
reach high orders in the expansion, allowing one to go to 
large N in Eq. (C3) and, correspondingly, to large t — to- 
The expansion coefficients are expressed in terms of 
Bessel functions of the first kind of order k. For more 
details, the reader is referred to Refs. [42, 43]. 

We use this approach to evaluate '^izi,Z 2 ,t) for t > 
—tr- We choose time steps around 0.2a;“^ and N up to 
800. A time step of 0.2a;“^ is small enough to resolve the 
tunneling dynamics and to extract the time dependence 
of the flux reliably, i.e., at the few percent accuracy level. 

Appendix D: Preparation of the initial state 

In Secs. Ill and IV we need the initial (equilibrium) 
state of the trapped particles at t = —R- We use 
p{t = —tr) = 0.795 for all cases. For this trap depth, 
tunneling is highly suppressed and the system is in a 
metastable state, i.e., it has a lifetime much larger than 
the time scale of the forthcoming tunneling process. To 
prepare the initial state, we “artificially” put a hard 
wall at z = lluho (the top of the barrier is located at 
z « Ooho)- We changed the position of the hard wall 
somewhat without seeing a notable change in the re¬ 
sults. For example, for the upper branch calculations 
at B = 900G, we changed the position of the hard wall 
to lOuho and 12aho and found that the overlap between 


the resulting initial states and the state prepared with 
the hardwall at lloho deviated from 1 by less than 10“®. 
This artificial boundary condition leaves the trap in the 
“inside” region unchanged and completely turns off the 
tunneling. The resulting eigenstates are to a very good 
approximation equal to the metastable states of the trap 
with finite barrier. 

We start with an initial wave packet that has a finite 
overlap with the state that we are looking for and act 
with the time-evolution operator, Eq. (C2) with imag¬ 
inary time T, on the initial wave packet [44, 45]. To 
propagate the wave packet in imaginary time, we use 
the real time-propagation methods discussed in Appen¬ 
dices C and E with t replaced by r/i. The initial wave 
packet can, in principle, be expanded in terms of un¬ 
known eigenfunctions of the Hamiltonian. After applica¬ 
tion of the time-evolution operator with imaginary time, 
each term in the expansion gets damped at a rate that is 
proportional to its energy. Thus, states with high energy 
decay fastest and eventually only the lowest energy state 
survives. We perform the imaginary time propagation 
using small At and normalize the wave packet to one 
after each step. This process can be generalized to find 
excited states by removing the lower energy eigenstates 
from the Hilbert space [45]. In practice, this is done by 
orthogonalizing the evolving wave packet and the lower 
energy eigenstate(s) after each time step. 

To implement the Chebychev expansion based ap¬ 
proach with imaginary time, we expand the exponential 
function in terms of real Chebychev polynomials and use 
the corresponding recursion relation [45]. We typically 
use about 15 terms in the series and time steps around 
(0.005w“^)/i. The Trotter formula based propagation 
scheme with imaginary time does not involve integrals 
over highly oscillatory functions and the calculations are 
computationally much less expensive than those for the 
real time propagation. 


Appendix E: Time propagation for Hamiltonian 
with two-body zero-range interaction 

The time propagation based on the Chebychev expan¬ 
sion is not applicable to the two-particle Hamiltonian 
with two-body zero-range interaction. In this case, we 
use a propagator that accounts for the two-body zero- 
range interaction exactly to determine the time evolution 
of the wave packet [46, 47]. This propagator has recently 
been used in Monte Carlo simulations for systems with 
zero-range interactions [48]. This appendix summarizes 
our implementation of the real time evolution in the pres¬ 
ence of a zero-range two-body potential. The wave packet 
41 at time t + At can be written as 

't/{zi,Z 2 ,t + ^t) = 

/ OO pOO 

/ pizi,Z2,zi,Z2;At)-^iz[,z!2,t)dz[dz2, (El) 

-OO J —OO 
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where the zero-range propagator p is defined through 
p{z[, z' 2 , Zi, Z2] l^t) = 

(zj, 2:2! exp(—02)- (E 2 ) 


In free space, i.e., when the two-body Hamiltonian con¬ 
sists of the kinetic energy and the zero-range interaction, 
the propagator pfree can be written as 


Pfree(4> 4; Zl,Z2\ At) = (z^ , Zi, At) 

XPfree(4, Z2, M)p'^l^{z[ - z'^, Zi - Z 2 , At), 


where 


1/2 


27TiAthJ 


exp 


p?L(z',z,At) = 
m(z — z')^' 
2iAth 


(E3) 


(E4) 


accounts for the single-particle kinetic energy and p^ 


,rel 


Pheeiz', 2 , At) = 1 - exp - 


m{zz' + | 2 :z'|) 
2iAth 


mm At piD 
4ti Y 


exp erfc(u), (E5) 


for the two-body zero-range potential [46, 47]. In 
Eq. (E5), erfc denotes the complementary error function 
and u is equal to m(|z| -I- \z'\+igiY)At/h)/V^miAth. For 
infinitely strong interaction, i.e., for |pid| = oo, Eq. (E5) 
simplifies to 


p[,t(z',z,At) 


l-exp(-^) for zz '>0 
0 for zz' < 0 


In the presence of the external potential we use the 
Trotter formula [49], 

p(z[,Z 2 ;zi,Z 2 ; At) « exp ^-^V;xt(4>4) 
xpfree( 4 . 4 ;^l)^ 2 ; At)exp (-^-^Ve^t(.Zi, Z 2 ) ) (E 6 ) 


This decomposition yields an error in the propagator that 
is proportional to At^. We use Eq. (El) with p given by 
Eq. (E 6 ) to propagate the wave packet in real time for 
each time step At. Unlike the Chebychev expansion ap¬ 
proach, the Trotter formula based approach is limited 
to small At. Importantly, the integrand in Eq. (El) os¬ 
cillates with a frequency that is proportional to I/At. 
To resolve these oscillations we need to choose a suffi¬ 
ciently dense spatial grid for the numerical integration 
of the right hand side of Eq. (El). We typically use a 
grid spacing Azj/oho < At/(10a;“^) (j = 1 and 2). We 
find that a value of At < 0.2a;“^ ensures that the norm 
of the wave packet, accounting for the absorbed portion 
of the wave packet, is 0.99 (or even closer to one) at the 
end of our simulation. Due to the need to evaluate the 
two-dimensional integral for each grid point, the Trotter 
formula based propagation scheme is much more com¬ 
putationally demanding than the Chebychev polynomial 
based propagation scheme. 


Appendix F: Application of the absorbing potential 


The damping of the wave packet in the numerical re¬ 
gions RiAni Ron, and RiBn has the same effect as an 
absorbing potential. After each time step, we multiply 
the wave packet by V{zi)V{z 2 ) [50], where 


V{z) 


1' 


f \ 

for 

z <Zd 

[ exp 

-a{ 

, A, / J 

for 

z> Zd 


Here, A^, a, and na are parameters whose values depend, 
in general, on the kinetic energy of the particle that is 
being absorbed. We use Ad = lOoho, a = 5, and Ud = 
2 with Zhw — Zd > 6 aho 7 where Zhw is the position of 
the hard wall at the end of the simulation grid. This 
parameter combination ensures that the reflection from 
the end of the numerical box is negligibly small. 


Appendix G: Flux analysis 


In this Appendix we discuss how to extract physical 
quantities from the density flux. P„(t) denotes the prob¬ 
ability to find n particles (n = 0 , 1 , or 2 ) inside the trap 
at time t. P 2 {t) is obtained by integrating the density 
|4'(zi, Z 2 , t)P over the region i ?2 (see Fig. 5), 


4 ^ 2 ( 1 )= [ |4'(zi,Z2,t)pfizidz2. (Gl) 

JR 2 

The initial condition is given by P 2 (—tr) = 1, i.e., at 
time t = —tr both particles are inside the trap. For 
t > —tr, we have P 2 it) + Pi{t) + Poit) = 1. The density 
\^{zi,Z 2 ,t)\'^ can flow from one region to another during 
the time propagation. To quantify the change of Pn{t), 
we use the current j(zi, Z 2 , t), 

j(zi,Z 2 ,t) = --^Im [4'*(zi,Z2,t) V4'(zi,Z2,t)] ,(G2) 

TO 

where 

_ 9 . 9 . 

V = —zi -f ;^z;2. (G3) 

UZ\ UZ2 


Here, Zi and Z 2 are the unit vectors in the Zi and Z 2 
directions, respectively. At each point in time and space, 
the continuity equation requires 


9|T(zi,Z2,t)P 

dt 


+ V • j(Zi,Z2,t) = 0. 


(G4) 


If we integrate Eq. (G4) over the region Ri {Ri can be 
equal to R 2 , Ria, or Rib', see Fig. 5), we find 


_9 



|4'(zi,Z2,t)P dzidz 2 = 



V • j(zi,Z 2 ,t) dzidz 2 . 


(G5) 


The left hand side of Eq. (G5) is the rate at which the 
probability of finding the system in region Ri changes. 
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To simplify the right hand side, we use the divergence 
theorem in two spatial dimensions, 


/ V ■ j{zi,Z2,t) dzidz2 = 
JRi 

]{zi,Z 2 ,t) ■ hi dl. 


(G6) 


Here, dl is the line element corresponding to the closed 
boundary Bi that encircles region Ri and is the unit 
vector perpendicular to the boundary and directed out 
of the region Ri. Equation (G 5 ) can thus be written as 


dt 



\'S{zi,Z 2 ,t)\'^ dzidz2 = 


-f jizi,Z2,t) - hi dl. 


Bi 


(G 7 ) 


The change of the probability to find the system in region 
Ri can be obtained from the flux through the boundary 
Bi. Applying Eq. (G 7 ) to region i?2, we obtain 


and 


1 AP2^l{t) 
P 2 {t) At 


(G 13 ) 


we have 72 = 7p + 7s- 7 p and 7 s oscillate in time for 
t not much larger than tr (typically t < 20ms) and are 
essentially constant for large t {t > 20 ms). The values 
reported in the main text are obtained by fitting the nu¬ 
merical data for sufficiently large t. 

If 2:d —>■ 00, we can hnd Pi{t) by integrating the density 
\'^{zi,Z2,t)\‘^ over the regions Ria and Rib (see Eig. 5 ) 
or, equivalently, by analyzing the flux through bound¬ 
aries b2AA , biA , 0 : ^2,IB, and 6ib,o- The average direc¬ 
tion of the flux is into the region Ria (-Rib) through 
boundary 62,1.4 (^2,ib) and out of the region Ria (Rib) 
through boundary biA,o (&ib,o)- In the upper branch 
simulations, we find vanishing flux through boundaries 
biA,o and 6ib,o- Thus, without worrying about the finite 
size of the simulation box, we can determine Pi (t) as the 
sum of the fluxes through boundaries 62,ia and 62,ib, 




' b2,lA 


j{zi,Z2,t) ■ hiA dl 


= ~ / j(^i>^2,t) • n2 d/ (G8) 
or 


R2(t) = 1 



j(zi, Z2, t) ■ n2 dl dt. 


B 2 


(G 9 ) 


To extract additional information from Eq. (G 7 ), we 
break the boundary B2 into pieces. In particular, flux 
through the boundary 62,0 corresponds to the corre¬ 
lated tunneling of two particles (pair tunneling) and flux 
through the boundaries 62,1^4 and 62,13 corresponds to 
single-particle tunneling (one particle tunnels and one 
remains in the trap). To quantify this in terms of tun¬ 
neling rates, we define the rate 72 at which P2it) decays 
during the time At through 


1 AP2{t) 
P 2 {t) At 


(GIO) 


Next, we divide the quantity AP2{t) into two pieces, 
namely the change AR2-s-o(0 due to the pair tunnel¬ 
ing (flux through the boundary 62,0) and the change 
AP2^i{t) due to the single-particle tunneling (flux 
through the boundaries 62,1,4 and 62,13), 


AP 2 {t) = AR2^o(i) + AP 2 ^i{t). (Gil) 


Defining the pair tunneling rate 7p and the single-particle 
tunneling rate 7s, 


1 AP2^o{t) 
P 2 {t) At 


(G 12 ) 


+ / i{zi,Z 2 ,t)-hiB dl\dt.{Gl 4 ) 

It should be noted that if the flux through boundaries 
61,4,0 and 613,0 is non-zero, then the evaluation of Pi{t) 
is more involved; this case is not discussed here. 


Appendix H: Additional comments on the WKB 
approximation 

As discussed in the main text, the WKB approx¬ 
imation yields single-particle tunneling rates that are 
smaller (larger) than the exact tunneling rates for the 
trap ground state (hrst excited trap state). To elabo¬ 
rate on this behavior, we determine p{t = 0) for the trap 
ground state, the first excited trap state, and the sec¬ 
ond excited trap state such that (a) T = 0.06267 and 
(b) T = 0 . 0063 . We then perform exact single-particle 
time propagation calculations for these cases, starting 
with a quasi-eigenstate (either the ground state, the first 
excited trap state, or the second excited trap state) for 
p{t = —tr) = 0 . 795 . Table VII summarizes the result¬ 
ing tunneling rates 7”p™. It can be seen that 7"p™ is 
approximately independent of the state number but, as 
expected, strongly dependent on the actual barrier the 
particle has to tunnel through. Due to the dependence of 
/WKB Qjj state (through the WKB energy), the WKB 
rates 7^^® for the three states vary by about a factor 
of 6 for cases (a) and (b). Eor the parameters considered 
in Table VB and in the main text, the WKB rate for 
the ground state is smaller than that obtained through 
the full time propagation, with the ratio 7^^^/7”p ™ de¬ 
pending on the exact shape of the trap. Eor the excited 
states, in contrast, the WKB rates are larger than those 
obtained through the full time propagation. 
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TABLE VII. Single-particle WKB versus exact tunneling rates. The tunneling coefficients for cases (a) and (b) are T = 0.06267 
and T = 0.0063, respectively. The third column reports the value of pit = 0) for which the trap ground state, first excited trap 
state, and second excited trap state yield the desired T. The fourth and fifth columns report the WKB frequency and the 

single particle WKB tunneling rate 7 ,^^®, Eq. ( 6 ), respectively. For comparison, the sixth column reports the tunneling rate 
7 “““ obtained from our exact time-dependent simulations. The calculations are performed for C = 1890G/m, Vo = 56.16i?ho, 
and zr = 8.548aho. 


case trap state p{t = 0 ) (ms 7 ^^® (ms 7 "™ (ms 


(a) 

gr. st 

0.63540 

0.322 

0.0202 

0.0330 

(a) 

1 st exc. 

st. 0.67687 

1.104 

0.0692 

0.0330 

(a) 

2 nd exc. 

st. 0.71486 

1.999 

0.1253 

0.0329 

(b) 

gr. st 

0.6489 

0.352 

0.00222 

0.00406 

(b) 

1 st exc. 

st. 0.6899 

1.1732 

0.00739 

0.00403 

(b) 

2 nd exc. 

st. 0.7277 

2.0965 

0.013215 

0.00400 
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